{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.1 Parabolic model problem\n",
    "We are solving the unsteady heat equation \n",
    "\n",
    "$$\\text{find } u:[0,T] \\to H_{0,D}^1 \\quad \\int_{\\Omega} \\partial_t u v + \\int_{\\Omega} \\nabla u \\nabla v + b \\cdot \\nabla u v = \\int f v  \\quad \\forall v \\in H_{0,D}^1, \\quad u(t=0) = u_0$$\n",
    "\n",
    "with a suitable advective field $b$ (the wind)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "import netgen.gui\n",
    "from math import pi\n",
    "from netgen.geom2d import SplineGeometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Geometry: $(-1,1)^2$\n",
    "* Dirichlet boundaries everywhere\n",
    "* Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle( (-1, -1), (1, 1), \n",
    "                 bcs = (\"bottom\", \"right\", \"top\", \"left\"))\n",
    "mesh = Mesh( geo.GenerateMesh(maxh=0.25))\n",
    "Draw(mesh)\n",
    "fes = H1(mesh, order=3, dirichlet=\"bottom|right|left|top\")\n",
    "\n",
    "u,v = fes.TnT()\n",
    "\n",
    "time = 0.0\n",
    "dt = 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We define the field $b$ (the wind) as \n",
    "\n",
    "$$b(x,y) := (2y(1-x^2),-2x(1-y^2)).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))\n",
    "Draw(b,mesh,\"wind\")\n",
    "from ngsolve.internal import visoptions\n",
    "visoptions.scalfunction = \"wind:0\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* bilinear forms for \n",
    " * convection-diffusion stiffness and \n",
    " * mass matrix seperately.\n",
    "* non-symmetric memory layout for the mass matrix so that a and m have the same sparsity pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes, symmetric=False)\n",
    "a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx\n",
    "a.Assemble()\n",
    "\n",
    "m = BilinearForm(fes, symmetric=False)\n",
    "m += u*v*dx\n",
    "m.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Implicit Euler method\n",
    "$$\n",
    "  \\underbrace{(M + \\Delta t A)}_{M^\\ast} u^{n+1} = M u^n + \\Delta tf^{n+1}\n",
    "$$\n",
    "\n",
    "or in an incremental form:\n",
    "\n",
    "$$\n",
    "  M^\\ast (u^{n+1} - u^n) = - \\Delta t A u^n + \\Delta tf^{n+1}.\n",
    "$$\n",
    "\n",
    "* Incremental form: $u^{n+1} - u^n$ has homogeneous boundary conditions (unless boundary conditions are time-dependent).\n",
    "* For the time stepping method: set up linear combinations of matrices.\n",
    "* (Only) if the sparsity pattern of the matrices agree we can copy the pattern and sum up the entries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "First, we create a matrix of the same size and sparsity pattern as m.mat:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": [
     "outputPrepend"
    ]
   },
   "outputs": [],
   "source": [
    "mstar = m.mat.CreateMatrix()\n",
    "print(m.mat.nze, a.mat.nze, mstar.nze)\n",
    "print(mstar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To access the entries we use the vector of nonzero-entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(mstar.nze)\n",
    "print(len(mstar.AsVector()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": [
     "outputPrepend"
    ]
   },
   "outputs": [],
   "source": [
    "print(mstar.AsVector())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the vector we can build the linear combination of the a and the m matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mstar = m.mat.CreateMatrix()\n",
    "mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()\n",
    "invmstar = mstar.Inverse(freedofs=fes.FreeDofs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We set the r.h.s. $f = exp(-6 ((x+\\frac12)^2+y^2)) - exp(-6 ((x-\\frac12)^2+y^2))$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = LinearForm(fes)\n",
    "gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))\n",
    "Draw(gaussp,mesh,\"f\")\n",
    "f += gaussp*v*dx\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and the initial data: $u_0 = (1-y^2)x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.Set((1-y*y)*x)\n",
    "Draw(gfu,mesh,\"u\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = gfu.vec.CreateVector()\n",
    "tstep = 1 # time that we want to step over within one block-run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "t_intermediate=0 # time counter within one block-run\n",
    "while t_intermediate < tstep - 0.5 * dt:\n",
    "    res.data = dt * f.vec - dt * a.mat * gfu.vec\n",
    "    gfu.vec.data += invmstar * res\n",
    "    t_intermediate += dt\n",
    "    print(\"\\r\",time+t_intermediate,end=\"\")\n",
    "    Redraw(blocking=False)\n",
    "print(\"\")\n",
    "time+=t_intermediate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Alternative version with iterative solvers\n",
    "\n",
    "* For a factorization of $M^\\ast$ (${M^\\ast}^{-1}$) we require a sparse matrix $M^\\ast$ \n",
    "* To store $M^\\ast$ as a sparse matrix requires new storage (and same memory layout of $A$ and $M$)\n",
    "* For iterative solvers we only require the matrix (and preconditioner) applications\n",
    "* `mstar = m.mat + dt * a.mat` has no storage but matrix-vector multiplications\n",
    "\n",
    "iterative solver version (with Jacobi preconditining):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mstar_alt = m.mat + dt * a.mat\n",
    "premstar_alt = m.mat.CreateSmoother() # + dt * a.mat.CreateSmoother()\n",
    "from ngsolve.krylovspace import CGSolver\n",
    "invmstar_alt = CGSolver(mstar_alt,premstar_alt, printrates='\\r',\\\n",
    "                        tol=1e-8, maxiter=200)\n",
    "\n",
    "print(premstar_alt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "tstep = 0.1\n",
    "t_intermediate=0 # time counter within one block-run\n",
    "while t_intermediate < tstep - 0.5 * dt:\n",
    "    res.data = dt * f.vec - dt * a.mat * gfu.vec\n",
    "    gfu.vec.data += invmstar_alt * res\n",
    "    t_intermediate += dt\n",
    "    # print(\"\\r t = {:24}, iteration steps: {}\".format(time+t_intermediate,invmstar_alt.GetSteps()), end=\"\")\n",
    "    Redraw(blocking=False)\n",
    "print(\"\")\n",
    "time+=t_intermediate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary material:\n",
    "* [Time dependent r.h.s.](#TDRHS)\n",
    "* [Time dependent boundary data](#TDBND)\n",
    "* [Runge-Kutta time integration](#RK)\n",
    "* [VTK Output](#VTK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 1: time-dependent r.h.s. data <a class=\"anchor\" id=\"TDRHS\"></a>\n",
    "\n",
    "Next: time-dependent r.h.s. data $f=f(t)$:\n",
    "\n",
    "* Use `Parameter` t representing the time. \n",
    "* A `Parameter` is a constant `CoefficientFunction` the value of which can be changed with the `Set`-function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = Parameter(0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "An example of a time-dependent coefficient that we want to use as r.h.s. in the following is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "omega=1\n",
    "gausspt = exp(-6*((x+sin(omega*t))*(x+sin(omega*t))+y*y))-exp(-6*((x-sin(omega*t))*(x-sin(omega*t))+y*y))\n",
    "Draw(gausspt,mesh,\"ft\",sd=2)\n",
    "time = 0.0\n",
    "from time import sleep\n",
    "while time < 10 - 0.5 * dt:\n",
    "    t.Set(time)\n",
    "    Redraw(blocking=False)\n",
    "    time += 1e-3\n",
    "    sleep(1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Accordingly we define a different linear form which then has to be assembled in every time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ft = LinearForm(fes)\n",
    "ft += gausspt*v*dx\n",
    "time = 0.0\n",
    "t.Set(0.0)\n",
    "gfu.Set((1-y*y)*x)\n",
    "#gfu.Set(CoefficientFunction(0))\n",
    "Draw(gfu,mesh,\"u\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "tstep = 10 # time that we want to step over within one block-run\n",
    "t_intermediate=0 # time counter within one block-run\n",
    "res = gfu.vec.CreateVector()\n",
    "while t_intermediate < tstep - 0.5 * dt:\n",
    "    t.Set(time+t_intermediate+dt)\n",
    "    ft.Assemble()\n",
    "    res.data = dt * ft.vec - dt * a.mat * gfu.vec\n",
    "    gfu.vec.data += invmstar * res\n",
    "    t_intermediate += dt\n",
    "    print(\"\\r\",time+t_intermediate,end=\"\")\n",
    "    Redraw(blocking=False)\n",
    "print(\"\")\n",
    "time+=t_intermediate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 2: Time dependent boundary conditions <a class=\"anchor\" id=\"TDBND\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* $u|_{\\partial \\Omega} = u_D(t)$, $f=0$\n",
    "* implicit Euler time stepping method, non-incremental form:\n",
    "\n",
    "  $$\n",
    "    M^\\ast u^{n+1} = (M + \\Delta t A) u^{n+1} = M u^n\n",
    "  $$  \n",
    "  \n",
    "* Homogenize w.r.t. to boundary conditions, i.e. we split \n",
    "\n",
    "  $$ u^{n+1} = u^{n+1}_0 + u^{n+1}_D $$ \n",
    "  \n",
    "  where $u^{n+1}_D$ is a (discrete) function with correct boundary condition:  \n",
    "  \n",
    "$$\n",
    "  {M^\\ast} u^{n+1}_0 = M u^n - {M^\\ast} u^{n+1}_D\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "uD = CoefficientFunction( [(1-x*x)*IfPos(sin(0.3*pi*t),sin(0.3*pi*t),0),0,0,0])\n",
    "time = 0.0\n",
    "t.Set(0.0)\n",
    "gfu.Set(uD,BND)\n",
    "Draw(gfu,mesh,\"u\")\n",
    "# visualization stuff\n",
    "from ngsolve.internal import *\n",
    "visoptions.mminval = 0.0\n",
    "visoptions.mmaxval = 0.2\n",
    "visoptions.deformation = 0\n",
    "visoptions.autoscale = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "tstep = 2 # time that we want to step over within one block-run\n",
    "t_intermediate=0 # time counter within one block-run\n",
    "res = gfu.vec.CreateVector()\n",
    "while t_intermediate < tstep - 0.5 * dt:\n",
    "    t.Set(time+t_intermediate+dt)\n",
    "    res.data = m.mat * gfu.vec\n",
    "    gfu.Set(uD,BND)\n",
    "    res.data -= mstar * gfu.vec\n",
    "    gfu.vec.data += invmstar * res\n",
    "    t_intermediate += dt\n",
    "    print(\"\\r\",time+t_intermediate,end=\"\")\n",
    "    Redraw(blocking=False)\n",
    "print(\"\")\n",
    "time+=t_intermediate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 3: Singly diagonally implicit Runge-Kutta methods <a class=\"anchor\" id=\"RK\"></a>\n",
    "\n",
    "\n",
    "We consider more sophisticated time integration methods, SDIRK methods. To simplify presentation we set $f=0$.\n",
    "\n",
    "SDIRK methods for the semi-discrete problem $\\frac{d}{dt} u = M^{-1} F(u) = -M^{-1} \\cdot A u$ are of the form:\n",
    "\n",
    "$$\n",
    "  u^{n+1} = u^n + \\Delta t M^{-1} \\sum_{i=0}^{s-1} b_i k_i\n",
    "$$\n",
    "\n",
    "with\n",
    "\n",
    "$$\n",
    "k_i = - A u_i \\quad \\text{ where $u_i$ is the solution to } \\quad \n",
    "(M + a_{ii} \\Delta t A) u_i = M u^n - \\Delta t \\sum_{i=0}^{i-1} a_{ij} k_j, \\quad i=0,..,s-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The coefficients a,b and c are stored in the butcher tableau:\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cccc}\n",
    "c_0 & a_{00} & 0 & \\ddots& \\\\\n",
    "c_1 & a_{10} & a_{11} & 0 & \\ddots \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & 0\\\\\n",
    "c_{s-1} & a_{s-1,0} & a_{s-1,1} & \\dots& a_{s-1,s-1} \\\\ \\hline\n",
    " & b_{0} & b_1 & \\dots& b_{s-1} \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "For an SDIRK method we have $a^\\ast = a_{ii},~i=0,..,s-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Simplest example: Implicit Euler\n",
    "$$\n",
    "\\begin{array}{c|c}\n",
    "1 & 1 \\\\ \\hline\n",
    " & 1  \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class sdirk1: #order 1 (implicit Euler)\n",
    "    stages = 1\n",
    "    a = [[1]]\n",
    "    b = [1]\n",
    "    c = [1]\n",
    "    astar = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can use for example the 2 stage SDIRK (order 3) method\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc}\n",
    "p & p & 0 \\\\\n",
    "1-p & 1-2p & p \\\\ \\hline\n",
    " & 1/2 & 1/2 \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "with $p = \\frac{3 - \\sqrt{3}}{6}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class sdirk2: #order 3\n",
    "    p = (3 - sqrt(3))/6\n",
    "    stages = 2\n",
    "    a = [[p, 0], \n",
    "         [1 - 2*p, p]]\n",
    "    b = [1/2, 1/2]\n",
    "    c = [p, 1 - p]\n",
    "    astar = p\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can use for example the 5 stage SDIRK (order 4) method\n",
    "$$\n",
    "\\begin{array}{c|cccc}\n",
    "1/4 & 1/4  \\\\\n",
    "3/4 & 1/2 & 1/4 & \\\\\n",
    "11/20 & 17/50 & -1/25 & 1/4 \\\\\n",
    "1/2 & 371/1360 & -137/2720 & 15/544 & 1/4 \\\\\n",
    "1 & 25/4 & -49/48 & 125/16 & -85/12 & 1/4 \\\\ \\hline\n",
    " & 25/4 & -49/48 & 125/16 & -85/12 & 1/4 \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class sdirk5: #order 4\n",
    "    stages = 5\n",
    "    a = [[1/4, 0, 0, 0, 0], \n",
    "         [1/2, 1/4, 0, 0, 0], \n",
    "         [17/50,-1/25, 1/4, 0, 0],\n",
    "         [371/1360, -137/2720, 15/544, 1/4,0],\n",
    "         [25/24, -49/48, 125/16, -85/12, 1/4]]\n",
    "    b = [25/24, -49/48, 125/16, -85/12, 1/4]\n",
    "    c = [1/4, 3/4, 11/20, 1/2, 1]\n",
    "    astar = 1/4    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### SDIRK2 :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "butchertab = sdirk2()   \n",
    "rhsi = gfu.vec.CreateVector()\n",
    "Mu0 = gfu.vec.CreateVector()\n",
    "ui = gfu.vec.CreateVector()\n",
    "k = [gfu.vec.CreateVector() for i in range(butchertab.stages)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We have to update the $M^\\ast$ matrix and reset initial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "time = 0.0\n",
    "t.Set(0.0)\n",
    "gfu.Set(uD,BND)\n",
    "Draw(gfu,mesh,\"u\")\n",
    "# visualization stuff\n",
    "from ngsolve.internal import *\n",
    "visoptions.mminval = 0.0\n",
    "visoptions.mmaxval = 0.2\n",
    "visoptions.deformation = 0\n",
    "visoptions.autoscale = 0\n",
    "\n",
    "mstar.AsVector().data = m.mat.AsVector() + butchertab.astar * dt * a.mat.AsVector()\n",
    "invmstar = mstar.Inverse(freedofs=fes.FreeDofs())\n",
    "invmass = m.mat.Inverse(freedofs=fes.FreeDofs())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "tstep = 2 # time that we want to step over within one block-run\n",
    "t_intermediate=0 # time counter within one block-run\n",
    "while t_intermediate < tstep - 0.5 * dt:\n",
    "    Mu0.data = m.mat * gfu.vec\n",
    "    for i in range(butchertab.stages):\n",
    "        # add up the ks as prescribed in the butcher tableau\n",
    "        rhsi.data = Mu0\n",
    "        for j in range(0,i):\n",
    "            rhsi.data += dt * butchertab.a[i][j] * k[j]\n",
    "        # Solve for ui (with homogenization for the boundary data)\n",
    "        t.Set(time+t_intermediate+butchertab.c[i]*dt)\n",
    "        gfu.Set(uD,BND)\n",
    "        ui.data = gfu.vec; rhsi.data -= mstar * ui\n",
    "        ui.data += invmstar * rhsi\n",
    "        # compute k[i] from ui\n",
    "        k[i].data = - a.mat * ui\n",
    "    t_intermediate += dt; t.Set(time+t_intermediate)\n",
    "    # Adding up the ks:\n",
    "    gfu.Set(uD,BND)\n",
    "    Mu0.data -= m.mat * gfu.vec\n",
    "    for i in range(0,butchertab.stages):\n",
    "        Mu0.data += dt * butchertab.b[i] * k[i]\n",
    "    gfu.vec.data += invmass * Mu0        \n",
    "    print(\"\\r\",time+t_intermediate,end=\"\")\n",
    "    Redraw(blocking=True)\n",
    "print(\"\"); time+=t_intermediate            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 4: VTK Output (\"exporting the nice pictures\") <a class=\"anchor\" id=\"VTK\"></a>\n",
    "\n",
    "* see also https://ngsolve.org/blog/ngsolve/2-vtk-output\n",
    "* or `py_tutorials/vtkout.py` (ngsolve repository)\n",
    "\n",
    "Outputting the nice pictures to vtk (to visualize with paraview):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "vtk = VTKOutput(mesh,coefs=[gfu],\n",
    "                names=[\"sol\"],\n",
    "                filename=\"vtk_example1\",\n",
    "                subdivision=3)\n",
    "vtk.Do()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "ls vtk_example1.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can also export vector fields:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "vtk = VTKOutput(mesh,coefs=[gfu,grad(gfu)],names=[\"sol\",\"gradsol\"],filename=\"vtk_example2\",subdivision=3)\n",
    "vtk.Do()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#%%bash\n",
    "#paraview vtk_example2.vtk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "And time dependent data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "vtk = VTKOutput(mesh,coefs=[gfu],names=[\"sol\"],filename=\"vtk_example3\",subdivision=3)\n",
    "gfu.Set((1-y*y)*x)\n",
    "vtk.Do(time=0)\n",
    "time = 0\n",
    "tstep = 1 # time that we want to step over within one block-run\n",
    "t_intermediate=0 # time counter within one block-run\n",
    "res = gfu.vec.CreateVector()\n",
    "i = 0\n",
    "while t_intermediate < tstep - 0.5 * dt:\n",
    "    res.data = dt * f.vec - dt * a.mat * gfu.vec\n",
    "    gfu.vec.data += invmstar * res\n",
    "    t_intermediate += dt\n",
    "    print(\"\\r\",time+t_intermediate,end=\"\")\n",
    "    Redraw(blocking=True)\n",
    "    i += 1\n",
    "    if (i%10 == 0):\n",
    "        vtk.Do(time=t_intermediate)\n",
    "print(\"\")"
   ]
  },
  {
   "source": [
    "Call paraview on gerenated data (if installed):"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "if ! command -v paraview &> /dev/null\n",
    "then\n",
    "    echo \"paraview is not installed.\"\n",
    "    exit\n",
    "else\n",
    "    paraview vtk_example3.pvd\n",
    "fi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.9.5 64-bit"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  },
  "interpreter": {
   "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}